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ABSTRACT 

We study the hypothesis that observed X-ray/extreme ultraviolet emission from 
coronae of magnetically active stars is entirely (or to a large part) due to the super- 
position of flares, using an analytic approach to determine the amplitude distribu- 
tion of flares in light curves. The flare-heating hypothesis is motivated by time series 
that show continuous variability suggesting the presence of a large number of super- 
imposed flares with similar rise and decay time scales. We rigorously relate the am- 
plitude distribution of stellar flares to the observed histograms of binned counts and 
photon waiting times, under the assumption that the flares occur at random and have 
similar shapes. Our main results are: (1) The characteristic function (Fourier trans- 
form of the probability density) of the expected counts in time bins At is <j>p(s,At) 
= expi-T^j^dtil - 4> a (sE(t,A))]}, where T is the mean flaring interval, <f> a ( s ) is 
the characteristic function of the flare amplitudes, and H(i, At) is the flare shape con- 
volved with the observational time bin. (2) The probability of finding n counts in time 
bins At is P c {n) = (2ir)~ 1 fQ 7T dse~' ins (f>F(s,At). (3) The probability density of photon 
waiting times x is Ps(x) = d x (f)F(i,x)/(r} with (r) = d x 4>F{i,x)\ x= o the mean count 
rate. An additive independent background is readily included. Applying these results 
to EUVE/DS observations of the flaring star AD Leo, we find that the flare amplitude 
distribution can be represented by a truncated power law with a power law index of 2.3 
±0.1. Our analytical results agree with existing Monte Carlo results of Kashyap et al. 
(2002) and Giidel et al. (2003). The method is applicable to a wide range of further 
stochastically bursting astrophysical sources such as cataclysmic variables, Gamma Ray 
Burst substructures, X-ray binaries, and spatially resolved observations of solar flares. 

Subject headings: Methods: statistical — Stars: activity — stars: flare — stars: individual 
(AD Leo) 



1. INTRODUCTION 

Since the late 1930s the solar corona has been known to be much hotter than the photosphere, 
although the latter delineates the surface through which energy enters the corona. The corona is 
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thus not in a state of simple heat flow equilibrium, and its high temperature must be supplied 
with some form of mechanical, electrical or radiative work which is dissipated at coronal heights. 
Several heating mechanisms have been proposed, such as Alfven (Hollweg 1983; Heyvaerts &; Priest 
1983; Marsh & Tu 1997) or ion cyclotron (Axford & McKenzie 1996) waves, or electric currents 
(Heyvaerts & Priest 1984) and magnetic reconnection (Parker 1972) converting magnetic shear into 
jets, and, ultimately, into heat. Some of these processes operate in a continuous mode, whereas 
others invoke transients. Comprehensive reviews may be found, e.g., in Narain & Ulmschneider 
(1990), Ulmschneider, Rosner & Priest (1991), Benz (1995), Narain & Ulmschneider (1996), and 
Priest & Forbes (2000). 

With increasing sensitivity of the observations it was recognized that the 'quiet' solar corona 
contains in fact 'microflares' and 'nanoflares' down to the actual resolution limit. Clearly, these 
events are the signature of non-equilibrium processes, and therefore represent promising candidates 
for the primary heating mechanism of the solar corona (Krucker &; Benz 1998; Benz & Krucker 
1999). 

Stars at high magnetic activity levels show properties that are reminiscent of continuous flaring 
processes, such as persistent coronal temperatures of 5-30 MK, high coronal electron densities, or 
continuous non-thermal radio emission. See Giidel et al. (2003) for a summary of observational 
results. Foremostly, X-ray and extreme ultraviolet (EUV) time series of such stars reveal a high 
level of flux variability on time scales similar to flares (Butler et al. 1986; Ambruster, Sciortino, & 
Golub 1987; Collura, Pasquini, & Schmitt 1988; Kashyap & Drake 1999; Audard, Giidel, & Guinan 
1999; Audard et al. 2000; Osten & Brown 1999). Limited detector sensitivity could give the wrong 
impression that the emission is due to steadily emitting coronal sources, while new, highly sensitive 
observations in at least some cases show little indication of a base level emission that could be 
defined as being genuinely 'quiescent' (Audard et al. 2003; Giidel et al. 2002). A coronal heating 
mechanism that should successfully explain X-ray and EUV emission from magnetically active stars 
inevitably has to address the question of continuous variability. 

In the context of microflare heating one is interested in the frequency distribution of flare 
amplitudes. In particular, it is essential to have reliable estimators of the small but frequent 
flares, since these presumably carry a major part of the total energy. However, these small flares 
substantially overlap, so that their amplitudes are hard to determine from the noisy observations. 
One way out is to use flare identification algorithms that also include flare superposition (Audard 
et al. 2000). Another way out is to resort to a forward model, which constrains the flare amplitude 
distribution by model assumptions. Within a Monte Carlo framework, this approach has been 
followed by Kashyap et al. (2002) and in Giidel et al. (2003). 

In the present article we explore analytical forward methods. We describe a statistical model of 
superimposed flares which provides an exact relation between the (univariate) probability density 
of the flare amplitudes and the (univariate) probability densities of the observed binned count rates 
and photon waiting times. As the method involves the histograms of binned count rates and photon 
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waiting times only, it is insensitive to permutation of the time bins and to the presence of long data 
gaps. The price for this robustness is the loss of temporal information, which is replaced by the 
assumption that all flares have similar shape and occur at random. See Isliker (1996) and Kraev & 
Benz (2001) for a complementary approach focusing on the temporal autocorrelation of the light 
curves (and thus on the first two moments of the univariate count rate distribution). 

In what follows we basically work with the observables 'counts' and 'count rates', with emphasis 
on the effect of flare superposition and counting noise. The conversion into luminosity and flare 
energy is briefly discussed in Section 6. 

2. MODEL 

While the energy release and relaxation processes are not known in detail, there is observational 
support for the similarity of flares on a global scale - at least for 'simple' flares. Stellar flares of the 
type investigated in this paper can usually well be represented by a sharp rise and an exponential 
decay. Examples supporting this view are found in de Jager et al. (1989), Giidel et al. (2002), Giidel 
et al. (2004), and Audard et al. (2003). The observed durations of flares, as investigated in the 
solar case, are almost independent of the flare amplitude or radiated energy even if 3-5 observed 
orders of magnitude are included. (Pearce & Harrison 1988; Shimizu 1995; Feldman et al. 1997; 
Aschwanden et al. 2000, see also summary in Giidel et al. 2003). The decay times of flares studied 
below were found to scatter around a constant value in Giidel et al. (2003), independent of peak 
flux. Since in stellar observations 3 orders of magnitude of flare energy are usually sufficient to 
explain the observed emission, we shall assume that a constant decay time is a good approximation 
for our case as well. For completeness, the case of a weak dependence of the decay time scale r on 
radiated total energy E, r oc E - 25 as suggested by one stellar investigation, was studied in both 
Kashyap et al. (2002) and Giidel et al. (2003), with the result that the main conclusions on the 
importance of small flares were reinforced. Observations also suggest that individual flares occur 
at random and independently of each other (for the solar case, see Wheatland 2000; Wheatland 
& Craig 2003). This independence might suggest that flares trigger an exhaustive 'reset' where 
mechanical stress (and thus memory) is completely destroyed. We shall thus adopt the following 
model for the instantaneous count rate in a fixed energy band: 

r(t) = J2 a kt{t-t k ) + b(t) [cts- 1 ] (1) 
k 

where £(t) represents the (dimensionless) flare shape, a k [ct s _1 ] are the flare amplitudes, 
and b(t) [ct s" 1 ] is a background intensity. We shall agree that < £(i) < £(0) = 1, and define 
r = f^Jiitydt to be the flare duration. All flare amplitudes a k are positive and drawn from the 
same probability density P a {x). (From here on, x denotes a generic real positive variable.) The 
flare times t k occur at random with mean flaring interval T [s], and it is assumed that all a k , t k , 
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and the background are statistically independent of each other. The observed counts are supposed 
to be a Poisson process with time- varying intensity r(t). Such a non-uniform Poisson process is 
not unambiguously defined, and we use here the construction shown in Fig. 1: the photon arrival 
times are distributed like the ^-components of uniformly distributed points (t, y) in the ty-plane 
satisfying y < r(t). Note that the background is included in r(t). When binned in time bins of 




t 

Fig. 1. — Construction of the time-dependent Poisson process with intensity r(t) from the standard 
2D Poisson process: imagine the whole 2D plane uniformly covered with points of density p = (r) 
x Is. By definition, the probability of finding no point in an arbitrary domain of area A is e~ pA . 
Points below the graph r(t) (boldface) are again a Poisson process, and so is their projection onto 
the t axis. The resulting t-coordinates are identified with the photon arrival times belonging to 
r(t). 

duration At, the counts are Poisson distributed with parameter (=expected number of counts in 
At) 



/t+At 
dt'r(t') 



(2) 
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Alternatively, we may work with the photon event list. In this case we consider the time 
intervals between consecutive photon arrivals ('photon waiting times'). 

Dimensional analysis shows that the flare-related part of the signal (1) is governed by the 
single dimcnsionless parameter r/T. The conversion into discrete time bins involves a second 
dimensionless parameter, At/r (or, alternatively, At/T). If the peaks are well-separated (r/T <C 1) 
and counting noise is low, one can directly estimate the flare amplitudes from the observed peaks in 
the light curve. If flares overlap (r/T > 1) or counting noise is relevant, this is no longer possible. 
We resort then to the model (1) and ask for its predictions on the observed probability densities of 
binned counts P c (n) and photon waiting times P$(x), which are available with high accuracy under 
stationary conditions. The goal is thus to establish the relation between P a (x) on the one hand, 
and P c (n) and P$(x) on the other hand. 

3. SOLUTION 

The linearity of Equation (1), together with the statistical independence of its ingredients, 
suggests that a rigorous mapping from P a (x) to P c {n) and P$(x) should exist. This is in fact the 
case, and the individual steps of this mapping are discussed in Sections 3.1 - 3.4. See Table 1 for 
an overview on our notation. 



Table 1: Notation") 





flare shape 




T 


flare duration [s], r = J^°^t^(t) 




T 


mean flaring interval [s] 




At 


observational time bin [s] 




~(t,At) 


= i: +A dt'm 




probability density 


quantity x, n 


characteristic function 


Pa(x) 


flare amplitudes [ct s -1 ] 


M s ) 


P f (x) 


flare intensity [ct s^ 1 ] 




Pf{x) 


flare parameter [ct] 


Ms, At) 


Pb(x) 


background intensity [ct s^ 1 ] 


M s ) 


Pb{x) 


background parameter [ct] 


<p B (s,At) 


P r (x) 


total intensity [ct 


Ms) 


Pr(x) 


total parameter [ct] 


Ms, At) 


P c (n) 


binned counts 


Ms, At) 


Ps(x) 


photon waiting time [s] 





a )The (Poisson-) 'intensity' is the instantaneous photon count rate; the (Poisson-) 'parameter' is the expected 
number of photons in time bins At. 'Total' means the sum of background and flares. 
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3.1. Flare Contributions 

We first ask for the probability that the flare signal f(t) = ^2 k a,k£,(t — tk) assumes values 
between x and x + dx. By statistical homogeneity, all times are equivalent. Consider thus the 
time t = 0, and let it be centered in an auxiliary interval / of length |/|. Neglecting flares outside 
I as well as the background, the intensity at t = is /(0) = ^2t k ei ak £( — Since /(0) is a 
sum of independent random variates, it is conveniently (Levy 1937) represented in a way which 
diagonalizes convolution, such as by its Laplace- or Fourier transform. We use here the latter, 
called the characteristic function (Lukacs 1970) and defined by 

POO 

Ms) = (e tsm )= / dxP f (x)e is * (3) 
J o 

where Pf(x) dx is the probability that /(0) is between x and x + dx. The probability density is ob- 
tained from the characteristic function by Fourier back transform, Pf(x) = J (fif(s)e~ tsx ds. 
We continue now with the calculation of the characteristic function <pf(s). The essential point to 
notice is that 4>f(s) can be partitioned according to the number of flares which occur in the interval 
I. One can thus write <f> f (s) = YlZo P l (e isfl ), where P, = e^l? {\I\/T) 1 /l\ is the probability of 
finding exactly I flares in /, and fi = Yl l k=i a k£,{—tk) is the corresponding photon intensity from 
these I superimposed flares in /. It is understood that /o = 0, because the case I = flare does 
not contribute to the count rate. Since the times tk are independent and uniformly distributed in 
/, the characteristic function of /; is 

(e lsfl ) = jf d ai P a ( ai ) ... d ai P a { ai ) jf\Jdh - Jf\l dt i exp{»*X>fc£(-t*)} ( 4 ) 

= J^ljdh ... Jdtttafsti-hj) ■...■<t> a (s£(-t l j) (5) 

= ^( fa </>*(* ■ (6) 

The step from equation (4) to equation (5) also invokes the fact that all ak are independently 
Pa-distributed, and the definition of the characteristic function of the flare amplitudes, 4> a {s) = 
(e zsa ). The step from equation (5) to equation (6) makes use of the symmetry of I about t = 0. 
Performing the sum over all possible numbers of flares and using equation (6), one finds that 



M*) = < e " /(0) > = E e "' J ' /r( | /|/T)t (e^') = exp { - T- 1 ! dt [1 - WW))]} ■ (7) 



We now want to liberate ourselves from the artificial time interval /, by extending it to infinity and 
thereby collecting all flare contributions. To this end we have to investigate whether the integral in 
equation (7) remains finite as |/| — > oo. This turns out to be the case under rather weak conditions, 
one of which is as follows (see Appendix A): if for some a > the flare shape £(t) decays faster 
than ji) -1 /" as \t\ — > oo, and if the probability density of the flare amplitudes decays no slower 
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than P a {x) ~ x~ l ~ a as x — > oo, then equation (7) remains well-behaved as |/| — > oo. Note that 
the above condition does not require the expectation value of the flare amplitudes to be finite; this 
freedom will allow us to consider every power law index v > 1 if the flare amplitudes have a power 
law distribution truncated at some lower cutoff, (v > 1 must be required to ensure that P a (x) can 
be normalized.) Letting thus / go to infinity, equation (7) becomes 

4> f (s) = exp { - T- 1 [°dt [1 - Ms C(t))]} • (8) 

J — oo 

The argumentation leading to equation (8) can be repeated, with minor modifications, for the 
Poisson parameter (expected counts) of time bins of duration At. Neglecting the background, the 
Poisson parameter of time bin [0, At] is F(0, At) = f Q t dt J2k a k£{t — tk) = J2k a k f-t* + 
so that we merely have to replace £(t) by 

pt+At 

~(t,At) = J dt'Z(t') (9) 
in order to obtain the characteristic function of F(0, At), 

^( S ,At) = ( e ^°' A *)) = exp{-T- 1 /'^ [l-^ a (sS(t,At))]}. (10) 

J — oo 

Equations (8,10) are our central result; they establish the connection between the probability 
density of flare amplitudes P a via <j) a and the distributions of photon intensities and expected 
counts, Pf, Pf via (pf, <pF- Equation (10) holds under the same conditions as equation (8). The 
function S(t, At) is a blurred version of the flare shape with resolution degraded to the observational 
time bin. It is easily seen that probabilistic normalization is preserved (0/(O)=l and 4>f{0, At)=l 
if 4> a (0)=l). It can also be shown that non-negative definiteness (Lukacs 1970) of 4> a {s) is inherited 
under the transformations (8) and (10), so that the Bochner (1959) theorem ensures that the 
Fourier inverses of (f>f(s) and ^>f(s, At) represent, in fact, probability densities. The transformation 
Pa( x ) ~^ Pf( x ) redistributes amplitudes x both towards higher (flare overlap) and lower (flare tails) 
values. It also broadens the probability density function by finite time binning, unless At C r, in 
which case E(t, At) -> At£(t) and <f> F (s,At) -> (f> f (sAt). In the limit At > T and if (a 2 ) < oo, 
normal behavior is approached, with F(0, At)/At concentrated about (a)r/T with relative spread 
of order (At/T) -1 / 2 . This expected result may be proven by formal expansion of equation (10) 
about 1/At = 0. 



3.2. Background 

So far the background has been neglected. When included, the total (flare plus background) 
photon Poisson intensity at time t = is r(0) = /(0) + 6(0), and the total Poisson parameter of the 
time bin [0, At] is R(0, At) = F(0, At) + 5(0, At) with B(0, At) = Jjffyt) dt. Since the background 
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is independent of the flares, the characteristic functions of the total photon Poisson intensity and 
-parameter are 

Ms) = Ms)4> f ( s ) (ii) 

ta(a,At) = M(s,At)M(s,At) (12) 

with (f>b{s) and cf>B(s,At) the characteristic functions of 6(0) and B(0,At), respectively. If the 
background is constant then 4>b(s, At) = (f) r (sAt) = e lsbAt . 



3.3. Binned Counts 

The counts in the time bin [0, At] are Poisson distributed with parameter R(0,At), where 
R(0, At) itself is PR-distributed. The binned counts have thus the characteristic function 

oo n —x ti r 

Ms, At) = (e isn ) = / dxP R (x) e -^-e isn = / dx P R {x) e^'^ = <j> R {i ~ ie is ,At) , (13) 
n=o J n. J 

which is 27r-periodic because the binned counts are integers. The probability density of the binned 
counts is obtained from Fourier back transform of 4> c (s, At); when expressed in terms of the original 
probability density P a (x), the probability of finding n counts in bins of duration At is 

P c (n) = —J M(s, At) exp {-ins- -J dt J dxP a (x)(l - e ~<^")^At)^ j (14) 

If a Gaussian approximation was used, with the binned counts modeled by c = R(0, At) + 
y/R(0, At) X with X standard normal, then c could assume any real number and the charac- 
teristic function is 4> c (s,At) ~ 4>r(s + is 2 /2,At). This approximation becomes exact at s — > 0, 
corresponding to large count rates. 



3.4. Photon Waiting Times 

Instead of time-binning the observed counts one may directly work with the photon event 
list, and consider the histogram of waiting times 5t between consecutive photon arrivals. This has 
the advantage that no ad-hoc choice of time bins is needed, and that the waiting time bins can 
be adapted for optimum (similar) probability content. In addition, long data gaps (3> dt) and 
detector saturation (St > T sa t) are easily segregated, and the construction of the photon waiting 
time histogram amounts to a direct determination of the moment generating function of r(t), which 
is (almost) the quantity provided by theory. The photon waiting times represent therefore a natural 
interface between theory and observation. 

To see why this is so, consider the 'waiting' event that there is (i) one photon between t\ 
and t\ + dt, (ii) no photon between t\ + dt and ti, and (iii) one photon between ti and ti + dt 
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(the waiting time is thus *2 — t\ > 0). By the binomial theorem, the probabilities of (i) and (hi) 
are r{t\)dt + 0(dt 2 ) and r(t2)dt + 0(dt 2 ), while from Figure 1 we see that the probability of 
(ii) is exp{— J)*V(t) dt}. Since (i), (ii) and (hi) are independent of each other, the waiting event 
occurs with probability density r(ti) exp{— J*/V(i) dt} rfo) = — c\<9i 2 exp{— / t *V(i) dt}. Denoting 
the average over realizations of r(t) by angular brackets, and using the translation invariance 
(exp{— / t *V(t) dt}} = (exp — / <2 ~*V(t) dt), one finds that the probability density of photon waiting 
times x = ti — t\ is 

P {W = J Vfl |( e -/.-K.)* ) = i_^^), (15 ) 

where the normalization constant N is determined by the constraint j^°Ps(x) dx = 1, implying 1 
that N = l/(r). While equation (15) is exact, it is generally difficult to calculate (j)ji(i,x) for 
arbitrary x. This is, however, not needed because the photon waiting time is much shorter than 
the flare duration in any practically relevant case. One may thus resort to the approximation 

4>n(i,x) ~ 4> r {ix) and 

W-"^ if and <a>r»l. (16) 

The cumulative waiting time distribution is, accordingly, 

r x Aj (ix) 

/ P 5 (y) dy ~ 1 + z-yy 2 if x < r and (°> r > 1 ■ ( 17 ) 
Jo \ r ) 

If the background is time- varying on time scale tj,, the conditions x <C Tfe and {b)rb 3> 1 must be 
added to the right hand side of equations (16, 17). The range of validity of equation (16) was 
explored analytically (Appendix B) and by Monte Carlo simulations (Section 4.2). It was found to 
hold even if the assumptions on the right hand side are not strictly valid. 

The result (16) can also be derived in the following, perhaps more intuitive, way: let P r (y) dy 
denote the probability that r{t) is between y and y + dy in a random time interval of fixed du- 
ration. Then, yP r (y)/ J y'P r (y')dy' is the probability that a given photon has been emitted at 
Poisson parameter y. (This probability is proportional to the total area of strips in the (t, y) 
plane with y < r(t) < y + dy - see Fig. 1). Assume now that r(t) varies only weakly during 
the local waiting time l/r(t). Thus yP r (y)/ f y'P r (y')dy' is also the probability that a wait- 
ing time 'belongs' to intensity y. Given the intensity y, the photon waiting time has proba- 
bility density ye~ y6t , so that the (marginal) probability density of photon waiting times x is 
P$(x) = f dy P r (y)y 2 e~ yx / f dyP r (y)y = —<f)"(ix)/(r), which is equation (16). 



"'"Because of J^°Ps(x) dx — N d x (exp— J^r(t) = — N lim^o d x (exp — f^rty) dt) — N{r). 
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3.5. Special Results 

For special flare shapes and amplitude distributions, equations (8) - (17) can be evaluated 
rigorously. This section summarizes results which are relevant for application to observational data 
(Section 5 below). 

3.5.1. Flare Shape 

Exponential decay is a generic feature of many relaxation processes, and is a natural model of 
flares. Observations of late-type main sequence stars often reveal a rather abrupt onset followed 
by an approximately exponential decay (for stellar examples, see Giidel et al. 2002, 2003; Audard 
et al. 2003), with decay time almost independent of flare size (e.g., Shimizu 1995; see also Section 
II). We shall thus use an exponential flare shape with constant decay time r as a proxy for exact 
calculations. If the flare shape is £(i) = 6(t) e~ l l T with 6(t) the Heaviside step function, then 

W = e*p{-l^I^M} (18) 

(see Appendix B). The abrupt and artificial onset of the one-sided exponential flare shape may 
be softened by an exponential increase, thus creating a two-sided exponential shape: £(t) = 
\9{— £)e*/ r+ + \Q{t)e~ t l T . Since flare contributions from t > and t < to £(0) are statisti- 
cally independent of each other, the characteristic function of their combined flare photon intensity 
is 

<f> f (s) = <l>J(s)c/>+( S ), (19) 

where 4>^(s) is given by equation (18) with r replaced by r^. Equations (18) and (19) show that 
two-sided exponential flares have the same characteristic function ^ as one-sided exponential flares 
with effective decay time r _ +r + . A similar result holds for <j>p if At <C r. From the observational 
point of view, this implies that one- and two-sided exponential flares cannot be distinguished based 
on their short-term photon waiting time (or binned count) distribution. We therefore assume in 
the sequel that the flare shape is one-sided exponential, with the understanding that our results on 
P a (x) are robust against the generalization of a one-sided to two-sided exponential shapes. 

3.5.2. Flare Amplitude Distribution 

For several forms of P a (x) the characteristic function (j)f( s ) can be explicitly obtained. Solvable 
cases include exponentials, power laws, and Levy distributions. Details are given in Appendix B; 
here we cite the results for power laws, which are a traditional and observationally well confirmed 
model of the distributions of solar (Datlowe, Elcan, & Hudson 1974; Lin et al. 1984; Krucker <fe 
Benz 1998; Benz & Krucker 1999; Aschwanden et al. 2000) and stellar (Collura et al. 1988; Audard 
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et al. 1999, 2000) flares. If P a (x) is a power law with exponent v, truncated at a lower cutoff x = ciq, 

P a (x) = (u-l)a u - 1 6(x-a )x~ u , (20) 

and if the flare shape is exponential, then (pf(s) can be exactly expressed in terms of special functions 
(Appendix B). For practical purposes, these may be approximated by elementary functions. Of 
particular interest is the cumulative photon waiting time distribution (eq. 17) which predicts the 
probability content of arbitrary photon waiting time bins: 

exp | - bx + \Ta<,x(c{a xy-l - ^ + 



7 . AraoQ-l) 
U ^ u-2 



(21) 



with C = tt\T(u) sin(7r^)] _1 . The relative accuracy of the approximation (21) was found to be better 
than 0.001 if a$x < 0.25 and v > 2.01. For v > 2 the expectation value of the flare amplitudes 
exists and is given by (a) = clq{v — — 2). In the limit v —* 2 the right hand side of equation 

(21) converges to the unit step function 6(x); this is also true if the exact equations (B10, Bll) are 
used. The limit (a) — > oo therefore predicts, theoretically, arbitrarily small photon waiting times. 



4. VALIDATION OF RESULTS 

Before application to real data, the results of Section 3 have been verified by several indepen- 
dent methods. The introduction of the characteristic function in Section 3.1 may appear somewhat 
arbitrary, and the question arises whether a more direct approach exists. The answer is yes, but the 
calculation is tedious and difficult to make rigorous. We shall therefore only sketch it in Appendix 
C. Here, we discuss validation based on an entirely independent analytical argument (Section 4.1) 
and on numerical simulations (Section 4.2). 



4.1. The Distribution of Consecutive Peaks 

For one-sided exponential flares - and only for these - the signal J2k a k£,(t — tk) is a continuous- 
time Markov process whose stationary measure can be obtained from a renewal argument. This 
provides an independent verification of Equation (18), and has some interest on its own as it allows 
to relax the Poisson assumption on the flare times (see Sect. 6 below). The argument is as follows. 
Discarding the background, the amplitudes at consecutive peaks are related by 

r n = a n + r n _ie-W T (22) 

with r n = r{t n + 0) and rj n = t n — t„_i the flare waiting times. Now, in the stationary case, the 
probability densities of r n and r„_i must agree. Therefore their characteristic function (j) Tn must 
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satisfy 

poo 

<Pr n ( S ) = (e lsr -) = < e ^)<e^- ie - WT ) = <f> a (s) / P v (x) (ae"*^) , (23) 

JO 

where Prj(x) is the flare waiting time distribution, and we have used the definition of the characteris- 
tic function and the fact that a n , r n and rj n are statistically independent of each other. For uniformly 
distributed flares, the flare waiting times are exponentially distributed, P v (x) = T~ l e~ x / T . Upon 
setting z = e~ v ^ T , equation (23) simplifies thus to 

<Pr n (s) = ^Ms)s~ T/T f z- l ^' T ^ rn {z)dz. (24) 
Equation (24) can be solved for <p Tn {s) by taking derivatives with respect to s, yielding 

KM = JJT)^ ~ T s + T MS) — ■ (25) 
Together with the initial (normalization) condition </v„(0) = 1, equation (25) has the unique solution 

A n (s) = Ms) exp { - T - f dy \ . (26) 

1 1 Jo y > 

The exponential factor on the right hand side of equation (26) acts as a low-pass filter, which 
broadens P r „(x) with respect to P a (x) by an amount proportional to t/T. As it must be, isolated 
flares reproduce the true flare amplitudes, lim T / r ^ 4>r„(s) = 4>a(s)- There is an apparent similarity 
between equations (18) and (26), and the exact relation is as follows. By construction, the peaks 
r n are the sum of two independent random variates: the value of the continuous signal just before 
the jump, r(t n — 0), plus the flare amplitude a n . From equation (26) we see that r(t n — 0) has 
the characteristic function exp{— rT _1 J s (l — (p a (x))x~ 1 dx}. Since r(t n — 0) occurs with equal 
probability than any r(t), the above expression must agree with the result of equation (18). This 
is in fact the case. 



4.2. Simulations 

Equations (8, 10, 14, 15) have also been verified by Monte Carlo simulations. In these simula- 
tions, a realization of r(t) is generated, and the binned Poisson parameters R(t, At) = f* +A r(t') dt' 
are numerically calculated. The binned counts are then simulated using a conventional Poisson 
number generator, and their histogram is created. For the photon waiting times we proceed as 
indicated in Figure 1: a large number iV = YT of points (t,y) are uniformly distributed in the 
rectangle [0, T] x [0,Y], with T ~ 3 x 10 3 r and Y ~ max r(t). The t-coordinates of points with 
y < r(t) constitute our photon list. Note that the background is included in r(t). The simulations 
are compared to numerical evaluations of eqns. (14,15,16) and to exact results (Appendix B). 

Figure 2 (left column) shows the theoretical (gray) and simulated (black) probability densities 
which occur in the transformation from the flare amplitudes into the binned counts. The flare 
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amplitudes are uniformly distributed between 1 and 3 ct s , the flare profile is a one-sided expo- 
nential with decay time r = 8 s, the mean flaring interval is T = 10 s, and the simulation involves 
2 x 10 6 time bins of duration At = 3 s. These parameters represent marginally resolved, partially 
overlapping flares. Accordingly, Pf(x) (Fig. 2 line 2 left) is smoothed with respect to P a (x) (line 
1 left) but has also a peak at x=0 which is due to low-count sequences between different flares. 
The binned background B(0, At) is modeled as chi square distributed with 3 degrees of freedom 
(Fig. 2 line 3 left); this ad-hoc choice yields a particularly simple characteristic function and is for 
illustration purposes only. The conversion of the binned Poisson parameter (Fig. 2 line 3 right) into 
discrete counts renders the characteristic function 27r-periodic (Fig. 2 line 4 right). For benchmark- 
ing purposes, the theoretical probability densities are obtained from numerical Fourier inversion 
of the associated characteristic functions (Fig. 2 right column, solid line - real part, dashed line - 
imaginary part). 

The transformation from the flare amplitudes into the photon waiting times is illustrated in 
Figure 3. The flare amplitudes are exponentially distributed (mean amplitude a = 0.15 ct 
mean flaring interval T = 20 s) and the flare shape is one-sided exponential (r = 40 s). The 
exponential flare amplitudes are chosen because they admit simple analytical results for Pg(x) (see 
Appendix B), which can be compared to the simulation. A constant background b = 0.1 ct 
is present, and the simulated photon list includes JV ~ 3 x 10 5 photons. The theoretical Pf(x) 
(Fig. 3 middle left gray) is a chi square distribution with r/T degrees of freedom (eq. B6), and 
Ps(x) (Fig. 3 bottom left gray) is given by the exact equation (B8). The characteristic functions of 
the flare amplitudes, flare photon intensity, and total (flare plus background) photon intensity are 
shown in the right column of Figure 3. Solid line denotes the real part, and dashed line denotes the 
imaginary part. A broad probability density corresponds to narrow characteristic functions, and 
vice versa. The constant background amounts to a modulation by e tsb . 

Both binned counts and photon waiting times have also been simulated for power law flare 
amplitudes, which are used as a fitting template for the EUVE/DS observations below. Figure 4 
shows the result of a simulation of one-sided exponential flares (r = 3000 s, T = 333 s) with lower 
cutoff ao = 0.007 ct and power law index v = 2.36. The simulation involves 3xl0 5 photons 
in 10 4 time bins of duration At = 100 s. Again, the black line represents simulation, and the gray 
line represents theory, with P$(x) being the short-term form (eq. 16) computed by differentiation 
of the exact expression (B10) for <fif(s). A constant background b = 0.09 ct s _1 is present, and an 
upper cutoff a± = 7 ct s" 1 has been applied in the simulations and when evaluating P c {n) from 
equation (14). 

5. OBSERVATIONS AND NUMERICAL RESULTS 

We have applied our theoretical results to a 3700 ks (elapsed time) observation of the active 
star AD Leo acquired by the Extreme Ultraviolet Explorer (EUVE, e.g. Malina & Bowyer 1991) 
Deep Survey (DS) instrument between April 2 and May 15 1999. This time series has previously 
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been analyzed by Kashyap et al. (2002) and in Giidel et al. (2003) by Monte-Carlo simulations of 
the model (eq. 1), so that we can compare our findings with existing results. The Deep Survey 
instrument records photons between 52 A and 246 A with maximum sensitivity at 91 A. Figure 5 
shows an overview of the light curve, time-integrated over spacecraft orbits (5663 s each). During 
one spacecraft orbit, AD Leo was in the EUVE field of view for typically 1500 s, followed by a 
~ 4200 s interval without observations. Furthermore, there are observation interrupts where the 
light curve of Figure 5 drops to zero. The period from 05/06 16:42 to 05/16 05:11 carries large 
contamination from energetic particles and is excluded from our analysis. 

Figure 6 shows the histograms of photon waiting times. The 'flares+background' are collected 
across a spatial window centered at AD Leo, while the 'background' is determined from a larger 
disjoint area. The resolution of our photon arrival times is 0.1 s. This is not the full EUVE/DS time 
resolution, but has been considered sufficient with regard to the actual count rates. Only photon 
waiting times from 0.25 s and 30 s are taken into account. The lower limit eliminates poor flare 
statistics and telemetry saturation at rare and extremely bright flares, and the upper limit eliminates 
low count rates and artifacts from (much longer) unobserved intervals. Our selection comprises 
118,637 photon waiting times, and the photon waiting time bin size is 0.1 s. The background waiting 
times (Fig. 6 'background' black line) are found to closely follow an exponential distribution (Fig. 
6 'background' gray line); accordingly, we model the background as constant. When scaled to the 
'flares+background' collecting area, the best-fit background intensity is (cf. Kashyap et al. 2002) 

b = 0.03 ± 0.002 ct s- 1 . (27) 

Once the background is known, we determine the flare amplitude distribution P a (x) under the 
assumption that P a {x) is a power law with exponent v and lower cutoff ao (eq. 20), and that the 
flare shape is one-sided exponential (Fig. 5 inlet). These assumptions correspond - up to a little 
relevant upper cutoff - to those made by Kashyap et al. (2002) and in Giidel et al. (2003). We also 
assert that the short-term forms (eqns. 16, 17) of the photon waiting time distribution apply. Our 
model parameters are thus (ao, v, t/T). They are determined by minimizing 

x 2 = E (ffn p Pn)2 ' ( 28 ) 

where H n (P n ) is the observed (predicted) number of waiting times in the n-th waiting time bin. 
The prediction is given by P n = N J* n+1 P$(x) dx/ Jq^ s Ps(x) dx, and we use the exact expression 
(B10) when computing the cumulative photon waiting time distribution (eq. 17). The sum in 
equation (28) runs over N = (30 — 0.25)/0.1 = 298 waiting time bins, and there are Adof = A Q — 
3 = 295 degrees of freedom. The minimum reduced chi square is 0.989, and the best-fit parameters 
are 

v = 2.29 ±0.07 a = 0.0049 ± 0.0015 ct s' 1 t/T = 11 ±4. (29) 

The resulting theoretical waiting time distribution is shown in Figure 6 ('flares+background' gray 
line) together with the observation ('flares+background' black line). The errors of (cio,v,t/T) are 
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not independent of each other, and the values given in equation (29) are projections of the 90% 
confidence domain (see Fig. 8 below). 

We also estimated the parameters (ao, v,t/T) from the observed histogram of binned counts 
(Fig. 7 black line). The (disjoint, adjacent) time bins have duration At = 100 s and contain, on 
average, 22 counts. Time bins which would intersect the transition from or to an unobserved time 
interval are rejected. In order to abate the arbitrariness of choosing a binning offset, the starting 
time of the first time bin has been varied at random within At, and the histogram of Figure 7 is 
an average over 50 random offsets. Time bins with less than 2 counts and more than 80 counts 
are excluded, because they are affected by the binning offset and have poor statistics. The best-fit 
solution is again found by minimizing the chi square expression of equation (28), where H n and 
P n are now the observed and theoretical frequencies of n counts in time bins of duration At. The 
theoretical prediction is P n = NP c (n) with TV = ^2 n H n = 4959 and n running over N Q = 80 — 2 = 
78 photon number bins. Since At <C r, we may use the short-term form of the probability density of 
binned counts, P c (n) = (2-k)- 1 ffis e- ins - bz ^ At <j) f (iz(s)At) with z(s) = l-e is and </> f (s) given by 
equation (B10). The minimum reduced chi square (per iVDOF = N Q — 3 = 75 degrees of freedom) 
is 1.003, which is attained at 

v = 2.31 ± 0.12 a = 0.0028 ± 0.0015 ct s' 1 t/T = 18 ± 7 . (30) 

The best-fit binned count distribution is superimposed in Figure 7 (gray line). The errors given in 
equation (30) are again projections of the 90% confidence domain. 

The confidence domains of {clq,v,t /T) for both the photon waiting time and binned count 
methods are visualized in Figure 8. The contours represent chi square sections through the best 
fit solution (crosses) at levels Xmin + (1-63, 2.33, 4.66) v / A r DOF, corresponding to (0.75, 0.9, 0.999) 
confidence levels 2 . Panels above the diagonal refer to the binned count method; panels below the 
diagonal refer to the photon waiting time method. The formal 90% confidence level is marked 
boldface. We find that the relative error of v is smaller than the errors of ao and t/T for both 
the binned count and waiting time methods. The power law index v can thus be estimated more 
reliably than the lower cutoff. The shape of the confidence domains indicates that the parameters 
are not independent of each other. The strongest interference is between the lower cutoff oo and 
the separation of individual flares, t/T, which reflects the fact that the average flare photon flux, 
tT~ 1 (cl) = rT~ 1 ao(l — v)/(2 — v), is sensitive to the product TT _1 ao only. Likewise, v is positively 
coupled to both ao and t/T because the function {y — \)/{y — 2) decreases with increasing v [y > 2). 
The results of the photon waiting time (eq. 29) and binned count (eq. 30) methods agree within the 
90% errors, which is considered as satisfactory, although the binned count method systematically 
yields lower ao and larger t/T than the photon waiting time method, whereas the product TT _1 ao 
is similar in both cases (0.0541 for the waiting time method; 0.0527 for the binned count method). 



2 If the model (ao, v, t/T) was true then the probability that the observed x 2 (eq. 28) deviates from its expectation 
value 7Vdof by more than e is about 1 — erf (e/2y/Nuo¥)- The cited 'confidence levels' represent the probabilities 
that this does not happen. 
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The mean flaring interval T and flare duration r enter the short-term forms of binned count- 
and photon waiting time distributions only in the dimensionless combination t/T. The decay 
constant t of flares with amplitudes > ao has been investigated by Kashyap et al. (2002) and in 
Giidel et al. (2003) and found to be about 3000 s; adopting this value, we find a mean flaring 
interval of 

T ~ 220 s (31) 

(average over both methods). The use of the short-term forms of the binned count- and photon 
waiting time distributions is justified, a posteriori, since At <C r, and since 5t < 30 s <C r and 
(a)r ~ 20 3> 1. In fact, direct numerical evaluations of the exact forms (eqns. 14, 15) at best-fit 
parameters (eqns. 29, 30, 31) deviate from their short-term approximations by less than 2%. 



6. DISCUSSION 

We have applied an analytical approach to the inversion of light curves to derive the amplitude 
distribution of constituent elementary flare light curves. The method uses either binned count 
light curves or the distribution of waiting times derived from photon lists. We have in particular 
considered flare shapes described by a one-sided exponential decay or a two-sided exponential 
profile, and we have studied flare amplitude distributions described by an exponential function and 
a truncated power law. The latter form has been favored in previous studies because solar flare 
statistics has revealed power-law distributions over a wide range in energy (Datlowe et al. 1974; Lin 
et al. 1984; Krucker & Benz 1998); explicit stellar studies indicate power laws as well (Collura et al. 
1988; Audard et al. 1999, 2000). We have thus used a power law to model EUVE/DS observations 
of AD Leo, where we find a power law index of v = 2.3 ± 0.1, a lower cutoff of ao = 0.004 ± 
0.002 ct and that there are t/T = 15 ± 8 flares per flare duration above ao; these values are 
compilations of the binned count- and photon waiting time results. 

Our results agree with complementary numerical methods recently applied to similar data sets 
(Audard et al. 1999, 2000; Kashyap et al. 2002; Giidel et al. 2003) indicating power-law indices 
for the flare amplitude distribution of typically 2.0 - 2.4, with most measurements concentrating 
around v = 2.2 - 2.3. While such values are necessary for power-law flare amplitude distributions to 
'explain' coronal heating on active stars, they are not sufficient. The emitted power observed from 
magnetically active coronae can only be due to the ensemble of flares if the power-law distribution 
extends to sufficiently low energies; this conjecture is of course difficult to prove due to the strong 
superposition of near-simultaneous flares that cannot be separated spatially, thus forming a pseudo- 
continuous emission level. When expressed in terms of flare energy, our numerical results (eqns. 
29, 30) correspond to a flare peak luminosity and minimum flare energy of 

a = (0.6...1.6) x 10 27 erg/s a r = (1.8 ... 4.8) x 10 29 erg , (32) 

where it is assumed that one EUVE/DS count corresponds to 2.8 x 10 29 erg emitted at AD Leo 
(see Giidel et al. 2003). Alternatively, a fraction of the observed emission may be truly steady (or 
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variable on time scales much longer than flare time scales), but this level should then be below 
the lowest level seen in light curves between individual detected flares; Audard, Giidel, & Skinner 
(2003) estimate that a mere 30% (or less) of the low-level emission of the dMe (binary) star UV 
Ceti could be ascribed to quasi-steady emission, the rest being recorded as variable on time scales 
of tens of minutes with shapes indicating continuous flaring. Additional information contained 
in the differential emission measure distribution could help interpreting the physical processes 
participating in the heating of the coronal plasma (Giidel et al. 2003). 

As mentioned in Giidel et al. (2003), a few words of caution are in order. A number of features 
cannot be recognized in the data per se and would require additional information. Specifically, 
the observations collect emission from across the visible hemisphere of the star, although heating 
mechanisms or, in our picture, flare amplitude distributions, may vary across various coronal struc- 
tures; further, some of the variability may be introduced by intrinsic changes in the large-scale 
magnetic field distributions; emerging or disappearing active regions will certainly contribute to 
variability, as do changes in the visibility of individual regions on the star as they move into or out 
of sight due to the star's rotation. Nevertheless, the light curve characteristics do not appear to be 
dominated by any of these longer-term effects, while variability on flare time scales is ubiquitous; 
similar conclusions were drawn by Audard et al. (2000), Kashyap et al. (2002), Giidel et al. (2003), 
and Audard et al. (2003). We therefore conclude that the present study further supports the view 
that stochastic flaring is of principal importance to coronae of magnetically active stars, and that 
coronae of such stars may significantly, if not entirely, be heated by flares. 

A practical benefit of the analytical method lies in its efficient assessment of the compatibility of 
model and observations. The computation of x 2 for a given set of parameters takes milliseconds on 
a medium-size work station, while a Monte Carlo simulation of comparable accuracy takes seconds 
or minutes due to slow (~ iV -1 / 2 ) convergence at rare large flares. The efficiency of the analytical 
approach allows to systematically explore the whole parameter space for the global chi square 
minimum. For instance, Figure 8 has been obtained by a brute- force calculation of x 2 ( a o,v, r/T) 
on a cubic lattice of dimension 90x95x100, taking about one hour computation time, whereas a 
similar Monte Carlo simulation would take many weeks. The possibility for systematic exploration 
of the parameter space (at low resolution) is important because the chi square function has generally 
several local minima, so that traditional minimization routines such as conjugate gradients may 
fail. Note that the simple contours of Figure 8 depict only the local behavior around the global 
minimum. 

When comparing the binned counts- and photon waiting time methods it was found that the 
binned counts method has larger confidence domains than the photon waiting time method, but is 
also more robust in finding the global minimum. A good strategy is thus to start with the binned 
counts, using a low-resolution exploration of the parameter space. As soon as a the approximate 
global chi square minimum is found, this may be refined by the photon waiting time method. In 
any case the final solutions should be cross-checked for both methods, to ensure consistency and to 
obtain an error estimate for the parameters which goes beyond the errors of the separate methods. 
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We would finally stress that our analytical method is in no way restricted to the particular 
time series studied in this paper. The method applies whenever the process can be described in the 
form of equation (1), which may include models for so divers phenomena as cataclysmic variables, 
Gamma Ray Burst substructures, and X Ray binaries. There is also no need to restrict the process 
(1) a the real line (see below). A drawback of the analytical method, compared to Monte Carlo 
simulations, is its limited capability to incorporate systematic instrumental artifacts. We close by 
mentioning a few straight forward extensions of the model (eq. 1). 

Flare clustering. The assumption of independent flare times might be questionable in the 
light of more advanced avalanche models (Vlahos et al. 1995) where flares may trigger each other. 
Flare clustering can be accounted for by a flare waiting time distribution which does not create 
uniform flare times {P^{x) 7^ T~ 1 e~ x / T ), but has a long-range tail permitting rare but large jumps. 
The 'clusters' are between the jumps. If the flare amplitudes are still independently P a -distributed, 
rise immediately and decay exponentially, then the functional Equation (23) for the peak values 
still holds, but the separation into the form (24) is no longer possible. However, Equation (23) can 
be iterated: the sequence 

4°]( S ) = Ms) 

POD 

4 k n +1) (s) = Ms) dxP v (x)4 k J(se-^) (33) 
J 

converges uniformly to the characteristic function of r(t n + 0) of clustered flares. All iterates are 
characteristic functions, and the characteristic function of the flare photon intensity is given by 
(j)f(s) = lim^oo^ (s)/M s ) ( see Section 4.1). 

CCD movies. The arguments leading to Equations (8) and (10) are not restricted to one- 
dimensional light curves. A useful generalization is to a sequence of spatially resolved CCD images 
od solar flares, where binning is over the space-time pixel A = Ax x At, and the auxiliary interval 
/ of Section 3.1 becomes I = L 2 x T. Provided that the flares occur independently and uniformly 
in /, the derivation of equation (10) goes through without modifications if At is replaced by A 
and E(t, At) is replaced by E(x, A) = f A £(x — x') dx' with x = (x, t) and £(x) the space-time flare 
shape. 

We thank Marc Audard for discussions on the EUVE data set. 



A. The limit |/| — > 00 

By definition, <j) a [s) satisfies 4> a (0) = 1 (normalization). Let us assume that 4> a (s) ~ 1 — 
c|s| a sign(s) as \s\ — ► with a > and c > 0; this is a valid form of (f) a (s) which is realized, 
for a < 2, in Levy distributions (Gnedenko & Kolmogorov 1954). Since £(i) — > as \t\ — ► 00, 
the integral in equation (7) remains well-behaved for all s when |/| —> 00 if J^° 00 {^(t)) a dt < 00, 
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requiring that £(t) decays faster than as \t\ — ► oo. So far about the decay of the flare shape. 

The decay of the probability density P a (x) at large flare amplitudes x is related by Tauberian 
theorems (Stein 1999) to the small-s limit of the characteristic function. For our assumption 
4> a (s) ~ 1 — c|s| a sign(s), one has that P a (x) = (2-7r)~ 1 ds e~ tsx (p a (s) ~ as x — ► oo. These 

are the conditions mentioned in Section 3.1. 



B. Explicit calculations 

For special flare shapes and amplitude distributions, equations (8) and (10) can explicitly be 
evaluated. To this end, it is advantageous to cast them in the Volterra form 



cf> f (s) = expj-T- 1 ^ f'dykfayHl-Mv))} ( B1 ) 

/•smaxH(t,A() 

b F (s,At) = expl-T- 1 ^/ dy Ki{ 8 ,y){l - <f> a (y))\ ( B2 ) 

7 is 



with kernels 



k{s,y)= Ki{s,y) = , (B3) 

s\^{ti)\ s|c/(tj,At)| 

where t\ is the Z-th solution of £(t) = y/s (kernel ki(s,y)) or E,(t, At) = y/s (kernel Ki(s,y)), 
respectively, and prime denotes derivative with respect to t. The success of an explicit evaluation 
of equations (Bl, B2) depends on the simplicity of 4> a (s) and on the ability to invert the functions 
£(i) and E(i, At) for the kernels ki(s,y) and Ki(s,y). (f>F(s,At) is generally harder to obtain than 

(f) f (s). 

If the flare shape is £(i) = #(t)e~' / '' r then the equation E(i, At) = y/s has the solutions t\ = 
—Tln(l — y/sT) — At (valid for the rise phase - At < t\ < 0) and t 2 = -rln(y/sr) + rln(l-e _A */ r ) 
(valid for the decay phase t 2 > 0). The corresponding kernel contributions (eq. B3) are K\ = 
t/(st — y) and K 2 = r/y, and equation (B2) becomes 

Ms,At) = e W -^ / dy(l-My))( + -)• (B4) 

T J \sT-y y) 

The mean of F(0, At) is T- 1 (a)rAt, and its variance is T^a 2 )^ (At - t - re" A '/ r ). For short 
time bins (At <C r) one has that y <C sr, and equation (B4) reduces to 4>F(s,At) = cf>f(sAt) 
with <j)f(s) given by equation (18). We present now several exactly solvable models for the flare 
amplitude distribution. 

For P a (x) = 9(x)a- 1 e- x / a and one-sided exponential flare shape, all characteristic functions 
and the photon waiting time distribution can be rigorously calculated. The results are 

<f> a (s) = (1-ias)- 1 (B5) 
c/) f ( s ) = (l-ias)- T / T (B6) 
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<t> F (s,x) = [{l-isaT)e x l T + isa?j ° Ta+l (B7) 

Equation (B6) indicates that Pf(x) is a chi square distribution with t/T (non- integer) degrees 
of freedom. The combination of exponential flare amplitudes and exponential flare profiles has 
the peculiar property that <f>f(s) = 4> a (s) (and thus Pf(x) = P a (x)) if t/T = 1. If a constant 
background b [ct s -1 ] is present, the exact (eq. 15) and short-term (eqns. 16, 17) forms of the 
photon waiting time distributions are 

t '"'((I • nr)i r or) n ' ,T ' 



-2 



Ps(x) 



b + arT- 1 
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2 n 2 Tp x/r ■ 
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e x/T -bar\ + ^ j (B8) 



„, , e- hx (l + ax)- T / T f/ arT' 1 ^ ^tT' 1 ^ 

Equations (B8) and (B9) allow to benchmark (and to be benchmarked by) our Monte Carlo sim- 
ulations (Section 4.2). They also provide a measure of validity of the short-term approximation 
(eqns. 16, 17): for example, for x = 0.33r, or = 25 photons per flare, t/T = 1 (moderate flare 
overlap) and b = 0.2arT _1 (20% background), the error of the short-term approximation is about 
0.6%. 

If P a (x) is a truncated power law (eq. 20) and if the flare shape is exponential, then 4>f{s) 
may be expressed in terms of special functions. Possible representations are 

ln^(a) = ^{ wT°T I + ia o s ^ 3J3QI, 1.2 - u], [2,2,3 - v],ia s)\ (BIO) 

M S) = e *(7-f+^)( aoS )-r/T e -^{E 1 (-iao S ) + K--o S )- 1 r(-,,- i a oS ) + ^} (BU) 

where p i^(n, d, ^) denotes the generalized hypergeometric function (Gradshteyn & Ryzhik 1980), 
7 = 0.5772 is Euler's constant, T(a,z) is the incomplete gamma function (Abramowitz & Stegun 
1970), Ei(z) = r(0, z) is the exponential integral, and the principal branches of the q F q and 
T functions are to be used. Equations (B10, Bll) have been obtained by the computer algebra 
program MAPLE. In numerical computations we use the form (B10) with a Taylor (small-argument) 
expansion of the hypergeometric function. For values of a$s <C 1, this expansion converges rapidly, 
and machine precision is reached after < 6 orders if aQS < 0.25. In fact, already order two yields 
< 0.0005 rms relative precision for Ps(x) over the whole waiting time range < 5t < 30s used in 
the EUVE data. This fact is exploited in the approximation (21), which is based on an order-two 
Taylor expansion of the hypergeometric function. 

If P a (x) is a Levy distribution (Levy 1937; Gnedenko Sz Kolmogorov 1954) with characteristic 
function 4> a (s) = exp{— c|s| a (l — ^jfj 2 ^)}, c > 0, and < a < 1, then P a (x) vanishes at negative 
x and therefore can model flare amplitudes. Large flare amplitudes (x » c) decay like P a {x) ~ 
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and superimposed flares have <fif(s) = (e 7 r/) _Ar / a |s|~ Ar exp{— ^Ei(ry|s| a )} with rj = c(l — 
ijlytan^) and 7 is Euler's constant (see above). Such Levy distributions might be a convenient 
model for small solar flares where power law indices < 2 have been reported (Crosby et al. 1993). 



C. Alternative Approach: Moments and Cumulants 

We ask for an elementary derivation of the results (8 - 15). The key observation is that 
equations (8) and (10) imply a simple relation between the moments of the flare amplitudes and 
the cumulants 3 of the instantaneous photon intensity and the binned Poisson parameter. This 
fact may be exploited to establish a relation between the moments of the flare amplitudes and 
the moments of the binned counts and photon waiting times, without recurrence to characteristic 
functions. We illustrate the procedure starting with the expectation value (/(0)). Interpreting the 
angular brackets as time average (which is justified by the assumption of stationarity) , the central 
limit theorem ensures that 

</(0)> = Bm^ij jTrft^Ofc^fc) = (a) ]^^Jdtm+0(Af} /2 ) - (a)± J £(t) & (CI) 

t k ei 

where A/j is the number of flares in the interval I (see Section 3.1). When restricting the first sum in 
eq. (CI) to flares inside I it has been supposed that contributions from outside I decay sufficiently 
fast to be neglected, and we have used the fact that all a& and t k are statistically independent of 
each other. Similarly, the second moment of /(0) is given by 



</(0) 2 > = lim ±- fdtfea&t-tkj) (2>£(*-<0) 

1 l_>0 ° 1 1 1 t k ei t t ei 

= { E + E i M , lim ^1 dt m (C2) 

= { E + E } m ^ jr\i dt i^fi jt\^ ~ tk) ^ ~ ti) 

- (a) 2 T- 2 [Jmdt) 2 + (a 2 )T~ 1 J[i{t)\ 2 dt. (C3) 

From equations (CI) and (C3) we recognize that (/(0) 2 ) - (/(0)) 2 = (a 2 )^ 1 f[£(t)} 2 dt: the 
variance of the instantaneous photon intensity equals the second moment of the flare amplitudes 
times a weight /[£(£)] 2 dt. This finding can be generalized to higher cumulants if the division into 



3 The cumulants of a quantity x are denned by n n (x) — (—i) n d™\n<t> x (s)\ s =o- Cumulants of order < n can be 
expressed in terms of moments of order < n, and vice versa; see, e.g., Eadie et al. (1971). 
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k=l and k^l contributions in equation (C2) is replaced by a sum over partitions into equal and 
non-equal indices. We shall not discuss here the (complicated) combinatorial weights, but note 
that they correspond to those of the cumulants expressed in terms of moments (e.g., Rota & Shen 
2000). As a consequence, 



Equation (C4) is a weak form of equation (8), since it represents its Taylor expansion around s=0. It 
is easily seen that equation (C4) also applies to the binned Poisson parameter F(0, At) if £(t) is re- 
placed by E(i, At). Therefore the cumulants (and thus the moments) of both /(0) and F(0, At) can 
be expressed in terms of the moments of the flare amplitudes. The probability density of the binned 
counts n is then given by P c (n) = (F{0, Ai) n e - F (°' Ai ))/n! = (n!)- 1 JX (-l) fc (F(0, Ai) n+fc )//fc!, 
while the probability density of short (Section 3.4) photon waiting times x is P$(x) oc {f(0) 2 e~ x ^ ^) = 
^2'kLo{—x) k (f(0) 2+k }/kl. Although these series may be evaluated in special cases, it is generally 
much more convenient to work with the characteristic functions, which collect all the moments in 
a compact way. 
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Fig. 2. — The transformation from a (rectangular) flare amplitude probability density (top left) to 
the probability density of the binned counts (bottom left). Black (gray) line indicate simulation 
(theory). The flare shape is exponential (t/T =0.8) and the bin size is At = 0.375r. Line 2 
(left) shows the probability density of the flare photon Poisson parameter; line 3 (left) shows 
the probability density of the binned background. The right column shows the real (solid) and 
imaginary (dashed) parts of the corresponding characteristic functions, except for line 3 (right) 
which refers to total (flare plus background) photon Poisson parameter. Arrows indicate the flow 
of the transformation P a (x) — > P c (n). 



-26- 




0.0 0.5 1.0 1.5 -40 -20 20 40 




0.0 0.5 1.0 1.5 2.0 -40 -20 20 40 




Fig. 3. — Similar to Figure 2, but for the photon waiting times. Left column: simulated (black) and 
theoretical (gray) probability densities of flare amplitudes (top), flare photon intensity (middle) 
and photon waiting times (bottom). Right column: characteristic functions of flare amplitudes 
(top), flare photon intensity (middle) and flare-plus-background photon intensity (bottom). The 
flare amplitudes are exponentially distributed with mean a=0.15 ct s _1 , the flare shape is one-sided 
exponential with r=40 s, the mean flaring interval is T=20 s, and there is a constant background of 
0.1 ct s -1 . The simulation involves 3xl0 5 photons. Arrows indicate the flow of the transformation 
P a (x) - P s {x). 
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Fig. 4. — Simulation: realization of r{t) (top left) of flares with one-sided exponential shape (r=3000 
s) and power law distributed amplitudes (top right), together with the histograms of binned counts 
and photon waiting times. Black (gray) line denotes simulation (theory). 
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Fig. 5. — EUVE/DS light curve from AD Leo obtained between April 2 and May 16 1999, integrated 
over one spacecraft orbit (5663 s). The model flare shape is shown as inlet, adopting the decay 
time of 3000 s of Kashyap et al. (2002). 
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Fig. 6. — Observed (black) and theoretical (gray) histograms of the photon waiting times of the 
EUVE/DS AD Leo observation (Fig. 5). The waiting time bin size is 0.1 s, the background is 
modeled as constant, and the flare amplitudes are modeled by a truncated power law. 
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Fig. 7. — Observed (black) and theoretical (gray) histograms of the binned counts of the EUVE /DS 
AD Leo observation of Figure 5. The time bin size is At = 100 s. 
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Fig. 8. — Chi square sections in parameter space (ao, v,t/T) through the best-fit solution (crosses). 
Panels above the diagonal refer to the binned count method; panels below the diagonal refer to the 
photon waiting time method. The contours represent (0.75, 0.9, 0.999) confidence levels. 



